/* Copyright (C) 2011 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

/***********************************************************
 * lammps.c -- functions for interfacing Brahms and LAMMPS * 
 ***********************************************************/

#include "dataStructs.h"

extern const VecR region, *siteCoords, *siteOrientVec; 
extern const Site *site;
extern const real timeNow, *mass, *charge, *dipole, *sigma, **sig, **eps;
extern const real springRigidity, angleRigidity, refAngle_CholPhosGly, refAngle_PhosGlyEst, refAngle_cisUnsat;
extern const int nSites, nLipids, nDOPCsDSPCs, nDOPEs, nTypes, nWaters;

// local variables:
static const int nAtomsPerLipid = 21; // assuming 15-site lipid models with 2 added atoms for each gly and est site
static const int nBondsPerLipid = 23; // assuming 15-site lipid models: 14 (main) + 6 (tip-tail charges) + 3 (orientation-restraining)

static const int nBondTypes = 9; // types defined below: 
static const int CHOL_PHOS = 1;
static const int PHOS_GLY = 2;
static const int GLY_EST =3;
static const int EST_CHAIN = 4;
static const int CHAIN_CHAIN = 5;
static const int AMINE_PHOS = 1;
static const int GLYPLUS_PHOS = 6;
static const int ESTPLUS_CHAIN = 7;
static const int GLYCHARGE_COM = 8;
static const int ESTCHARGE_COM = 9;

static const int nAnglesPerLipid = 16; // assuming 15-site lipid models: 13 (main) + 3 (atomic charges)
static const int nAngleTypes = 5; // 1_chol(amine)-phos-gly, 2_phos-gly-est, 3_PI_sat, 4_cisUnsat, 5_atomicCharges
static const int _phos_ = 1;
static const int _gly_ = 2;
static const int _piSat_ = 3;
static const int _cisUnsat_ = 4;
static const int _piDipCharges_ = 5;


/*******************************************************************************************************
 * WriteLammpsData -- writes "data." file to be read into the LAMMPS program through read_data command *
 *******************************************************************************************************/
void WriteLammpsData( FILE *fp )
{
  int n, id, atomCount, bondCount, angleCount, atom1, atom2, atom3;
  int idPhos, idGly, idEstA, idEstB, idGlyDipTip, idGlyDipTail, idEst, estDipTail, estDipTip, idChain1, idChainI, idChainJ;
  int nAtoms = nWaters + nAtomsPerLipid*nLipids; // #atoms for LAMMPS topology
  int nBonds = nBondsPerLipid*nLipids;
  int nAngles = nAnglesPerLipid*nLipids;
  int nAtomTypes;
  real disp, atomCharge;
  VecR tailCoords, tipCoords;

  if ( nLipids == nDOPCsDSPCs || nLipids == nDOPEs ) nAtomTypes = nTypes-1; // bilayer containing only one lipid species
  else nAtomTypes = nTypes; // mixed bilayer

  fprintf(fp, "LAMMPS 'data.' description - created with Brahms (www.personal.soton.ac.uk/orsi/brahms)\n");
  fprintf(fp, "\n");

  fprintf(fp, "%9d atoms\n", nAtoms);
  fprintf(fp, "%9d bonds\n", nBonds);
  fprintf(fp, "%9d angles\n", nAngles);
  fprintf(fp, "\n");
  fprintf(fp, "%9d atom types\n", nAtomTypes);
  fprintf(fp, "%9d bond types\n", nBondTypes);
  fprintf(fp, "%9d angle types\n", nAngleTypes);
  fprintf(fp, "\n");

  fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.x*LENGTH_Angstrom/2, region.x*LENGTH_Angstrom/2, "xlo xhi");
  fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.y*LENGTH_Angstrom/2, region.y*LENGTH_Angstrom/2, "ylo yhi");
  fprintf(fp, "% 13.7G% 13.7G%13s\n", -region.z*LENGTH_Angstrom/2, region.z*LENGTH_Angstrom/2, "zlo zhi");
  fprintf(fp, "\n");
  
  fprintf(fp, "Atoms\n\n");
  /*   DO_SITE { */
  /*     fprintf(fp, "%7d%5d", n+1, site[n].type+1 ); // write atom-ID and atom-type */
  /*     fprintf(fp, "%11.4f%11.4f%11.4f", siteCoords[n].x*LENGTH_Angstrom, siteCoords[n].y*LENGTH_Angstrom, siteCoords[n].z*LENGTH_Angstrom); */
  /*     fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID (this is 0 for water sites) */
  /*     fprintf(fp, "%7.2f", charge[site[n].type]); // write charge (this is 0.00 for water sites) */
  /*     fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole component ('dipole' style) */
  /* 	    siteOrientVec[ n ].x * dipole[site[n].type] * DIPOLE_D * D_IN_eA, */
  /* 	    siteOrientVec[ n ].y * dipole[site[n].type] * DIPOLE_D * D_IN_eA, */
  /* 	    siteOrientVec[ n ].z * dipole[site[n].type] * DIPOLE_D * D_IN_eA ); */
  /*     // write attributes for 'sphere' style - diameter & atom density: */
  /*     fprintf(fp, "%5.1f%5.1f", sigma[site[n].type]*LENGTH_Angstrom, mass[site[n].type] / (4*PI/3 * Cube((sigma[site[n].type]/2)*LENGTH_Angstrom))); */
  /*     fprintf(fp, "\n"); */
  /*   } */
  atomCount = 0;
  DO_SITE {
    atomCount++;
    if ( site[n].type==AMINE_TYPE )
      fprintf(fp, "%7d%5d", atomCount, 2 ); // write atom-ID and atom-type foe amine - TO SORT!
    else fprintf(fp, "%7d%5d", atomCount, site[n].type+1 ); // write atom-ID and atom-type
    fprintf(fp, "%11.4f%11.4f%11.4f", siteCoords[n].x*LENGTH_Angstrom, siteCoords[n].y*LENGTH_Angstrom, siteCoords[n].z*LENGTH_Angstrom); 
    fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID (this is 0 for water sites, 1 for lipid-1, 2 for lipid-2, etc.)
    fprintf(fp, "%7.2f", charge[site[n].type]); // write charge (this is 0.00 for point-dipole sites):
    if (site[n].type==GLYCEROL_TYPE || site[n].type==ESTER_TYPE) { // generate two extra charged atoms to represent the dipole
      fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set dipole to zero as it is going to be represented by two charges (below)
    }
    else { fprintf(fp, "% 10.5f% 10.5f% 10.5f", // write dipole components (non-zero only for water sites)
		   siteOrientVec[ n ].x * dipole[site[n].type] * DIPOLE_D * D_IN_eA, 
		   siteOrientVec[ n ].y * dipole[site[n].type] * DIPOLE_D * D_IN_eA, 
		   siteOrientVec[ n ].z * dipole[site[n].type] * DIPOLE_D * D_IN_eA );
    }
    // write attributes for 'sphere' style - diameter & atom density:
    fprintf(fp, "%5.1f%5.1f", sigma[site[n].type]*LENGTH_Angstrom, mass[site[n].type] / (4*PI/3 * Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
    fprintf(fp, "\n"); 
    if (site[n].type==GLYCEROL_TYPE || site[n].type==ESTER_TYPE) { // generate two extra charged atoms to represent the dipole
      disp = sigma[site[n].type] / 4; // displacement from dipole center
      VSAdd( tipCoords, siteCoords[n], disp, siteOrientVec[n] );
      VSAdd( tailCoords, siteCoords[n], -disp, siteOrientVec[n] );
      atomCharge = dipole[site[n].type] * DIPOLE_D * D_IN_eA / (2*disp*LENGTH_Angstrom);
      // add "tip" charge atom:
      atomCount++;
      fprintf(fp, "%7d%5d", atomCount, CHARGE_TYPE ); // write atom-ID and atom-type
      fprintf(fp, "%11.4f%11.4f%11.4f", tipCoords.x*LENGTH_Angstrom, tipCoords.y*LENGTH_Angstrom, tipCoords.z*LENGTH_Angstrom); 
      fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID 
      fprintf(fp, "%7.2f", atomCharge); // write charge
      fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set point-dipole to zero 
      fprintf(fp, "%5.1f%5.1f", sigma[site[n].type]*LENGTH_Angstrom, mass[site[n].type] / (4*PI/3 * Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
      fprintf(fp, "\n"); 
      // add "tail" charge atom:
      atomCount++;
      fprintf(fp, "%7d%5d", atomCount, CHARGE_TYPE ); // write atom-ID and atom-type
      fprintf(fp, "%11.4f%11.4f%11.4f", tailCoords.x*LENGTH_Angstrom, tailCoords.y*LENGTH_Angstrom, tailCoords.z*LENGTH_Angstrom); 
      fprintf(fp, "%7d", site[n].inLipid ); // write molecule-ID 
      fprintf(fp, "%7.2f", -atomCharge); // write charge
      fprintf(fp, "% 10.5f% 10.5f% 10.5f", 0.,0.,0.); // set dipole to zero 
      fprintf(fp, "%5.1f%5.1f", sigma[site[n].type]*LENGTH_Angstrom, mass[site[n].type] / (4*PI/3 * Cube((sigma[site[n].type]/2)*LENGTH_Angstrom)));
      fprintf(fp, "\n"); 
    }
  }
  if (atomCount!=nAtoms) {printf("Error: atomCount!=nAtoms.\n"); exit(0);}
  fprintf(fp, "\n"); 
  
  fprintf(fp, "Bonds\n\n");
  bondCount = 0;
  DO_SITE {
    id = n + (site[n].inLipid - 1) * 6; // conversion to LAMMPS atom-ID
    switch ( site[ n ].lipidUnit ) {
    case WATER:     /* break - nonlipid site */ 
    case TAIL_A5: // break - this case is considered in TAIL_A4
    case TAIL_B5: // break - this case is considered in TAIL_B4
      break;
    case CHOLINE_OR_AMINE: /* choline***phosphate (DOPC) or amine***phosphate (DOPE) */
      if (site[n].type == CHOLINE_TYPE) { /* choline***phosphate (DOPC) */
	atom1 = id + 1;
	atom2 = atom1 + 1;
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHOL_PHOS, atom1, atom2 );
      } else { /* amine***phosphate (DOPE) */
	atom1 = id + 1;
	atom2 = atom1 + 1;
	fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, AMINE_PHOS, atom1, atom2 ); 
      }
      break;
    case PHOSPHATE: // phosphate***glycerol 
      atom1 = id + 1;
      atom2 = atom1 + 1;
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, PHOS_GLY, atom1, atom2 ); 
      break;
    case GLYCEROL: 
      idGly = id + 1; // gly
      idEstA = idGly + 3; // because of additional 2 sites used for glycerol atom charges
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstA ); // GLYCEROL***ESTER_A:
      idEstB = idGly + 11; // because of additional 4 sites used for glycerol and ester atom charges
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLY_EST, idGly, idEstB ); // GLYCEROL***ESTER_B:
      idGlyDipTip = idGly + 1;
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLYCHARGE_COM, idGly, idGlyDipTip ); // Gly COM *** dipole tip ('+' end) 
      idGlyDipTail = idGly + 2;
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLYCHARGE_COM, idGly, idGlyDipTail ); // Gly COM *** dipole tail ('-' end)
      idPhos = id; // phosphate (assumes phosphate'ID always precedes glycerol's ID
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, GLYPLUS_PHOS, idPhos, idGlyDipTip ); // Gly dipole tip ('+' end) *** phosphate
      break;
    case ESTER_A: // ESTER_A***TAIL_A1
    case ESTER_B: // ESTER_B***TAIL_B1
      if (site[n].lipidUnit == ESTER_A) idEst = id + 3;
      else idEst = id + 5;
      idChain1 = idEst + 3; 
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, EST_CHAIN, idEst, idChain1 ); 
      estDipTail = idEst + 2; 
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, ESTPLUS_CHAIN, estDipTail, idChain1 ); 
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, ESTCHARGE_COM, idEst, estDipTail ); // Est COM *** dipole tail
      estDipTip = idEst + 1; // est tip
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, ESTCHARGE_COM, idEst, estDipTip ); // Est COM *** dipole tip
      break;
    case TAIL_A1: // TAIL_A1***TAIL_A2
    case TAIL_A2: // TAIL_A2***TAIL_A3
    case TAIL_A3: // TAIL_A3***TAIL_A4
    case TAIL_A4: // TAIL_A4***TAIL_A5
      idChainI = id + 5;
      idChainJ = idChainI + 1;
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI, idChainJ ); 
      break;
    case TAIL_B1: // TAIL_B1***TAIL_B2
    case TAIL_B2: // TAIL_B2***TAIL_B3
    case TAIL_B3: // TAIL_B3***TAIL_B4
    case TAIL_B4: // TAIL_B4***TAIL_B5
      idChainI = id + 7;
      idChainJ = idChainI + 1;
      fprintf(fp, "%7d%5d%8d%8d\n", ++bondCount, CHAIN_CHAIN, idChainI, idChainJ ); 
      break;
    default:
      printf("Error: default case reached in WriteLammpsData. Exiting...\n");
      exit(0);
    }    
  }
  if (bondCount!=nBonds) {printf("Error: bondCount=%d!=%d=nBonds.\n",bondCount,nBonds); exit(0);}
  fprintf(fp, "\n");

  fprintf(fp, "Angles\n\n");
  angleCount = 0;
  DO_SITE {
    id = n + (site[n].inLipid - 1) * 6; // conversion to LAMMPS atom-ID
    switch ( site[ n ].lipidUnit ) { /* identify the site "n" corresponding to the angle "(n-1)--(n)--(n+1)" */
    case WATER: 
    case CHOLINE_OR_AMINE:
    case TAIL_A5: 
    case TAIL_B5: 
      break; /* there is no angle potential centered on these sites */
    case PHOSPHATE: /* choline_or_amine(1)-phosphate(2)-glycerol(3) */
      atom1 = id; // choline or amine
      atom2 = id+1; // phosphate
      atom3 = id+2; // glycerol
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _phos_, atom1, atom2, atom3 ); 
      break;
    case GLYCEROL:
      atom1 = id; // phosphate
      atom2 = id+1; // glycerol
      atom3 = id+4; // ester_a
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2, atom3 ); // phos-gly-estA
      atom3 = id+12; // ester_b
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _gly_, atom1, atom2, atom3 ); // phos-gly-estB
      atom1 = id+2; // dipTail 
      atom2 = id+1; // glycerol
      atom3 = id+3; // dipTip
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piDipCharges_, atom1, atom2, atom3 ); /* dipTail-glyCOM-dipTip: */
      break;
    case ESTER_A: 
      atom1 = id; // gly
      atom2 = id+3; // ester_a
      atom3 = id+6; // tail_a1
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); /* glycerol(3)-ester_A(4)-hydrocarbon_a1(5): */ 
      atom1 = id+4; // dipTail 
      atom2 = id+3; // ester_a
      atom3 = id+5; // dipTip
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piDipCharges_, atom1, atom2, atom3 ); /* dipTail-estCOM-dipTip: */ 
      break;
    case TAIL_A1:   
      atom1 = id+2; // estA
      atom2 = id+5; // chain_a1
      atom3 = id+6; // chain_a2
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); /* ester_A(4)-hydrocarbon_a1(5)-hydrocarbon_a2(6) */
      break;
    case TAIL_B1:   /* ester_B(4)-hydrocarbon_b1(5)-hydrocarbon_b2(6) */
      atom1 = id+4; // estB
      atom2 = id+7; // chain_b1
      atom3 = id+8; // chain_b2
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); 
      break;
    case TAIL_A3:   /* hydrocarbon_a2(6)-hydrocarbon_a3(7)-hydrocarbon_a4(8) */
    case TAIL_A4:   /* hydrocarbon_a3(7)-hydrocarbon_a4(8)-hydrocarbon_a5(9) */
      atom1 = id+4; 
      atom2 = atom1+1;
      atom3 = atom2+1; 
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); 
      break;
    case TAIL_B3:   /* hydrocarbon_b2(6)-hydrocarbon_b3(7)-hydrocarbon_b4(8) */
    case TAIL_B4:   /* hydrocarbon_b3(7)-hydrocarbon_b4(8)-hydrocarbon_b5(9) */
      atom1 = id+6; 
      atom2 = atom1+1;
      atom3 = atom2+1; 
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); 
      break;
    case ESTER_B: /* glycerol(3)-ester_B(10)-hydrocarbon_b1(11) -- SPECIAL CASE */
      atom1 = id-6; 
      atom2 = id+5;
      atom3 = id+8; 
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piSat_, atom1, atom2, atom3 ); 
      atom1 = id+6; // dipTail 
      atom2 = id+5; // ester_b
      atom3 = id+7; // dipTip
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _piDipCharges_, atom1, atom2, atom3 ); /* dipTail-estCOM-dipTip: */ 
      break;
    case TAIL_A2:   /* hydrocarbon_a1(5)-hydrocarbon_a2(6)-hydrocarbon_a3(7) */
      atom1 = id+4; 
      atom2 = atom1+1;
      atom3 = atom2+1; 
      // set cisUnsat angle type for DOPC and DOPE - case for DSPC to be coded!
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _cisUnsat_, atom1, atom2, atom3 ); 
      break;
    case TAIL_B2:   /* hydrocarbon_b1(5)-hydrocarbon_b2(6)-hydrocarbon_b3(7) */
      atom1 = id+6; 
      atom2 = atom1+1;
      atom3 = atom2+1; 
      // set cisUnsat angle type for DOPC and DOPE - case for DSPC to be coded!
      fprintf(fp, "%7d%5d%8d%8d%8d\n", ++angleCount, _cisUnsat_, atom1, atom2, atom3 ); 
      break;
    default:
      printf("Error in function WriteLammpsData - default case reached. Exiting...\n"); exit(0);
    } 
  }
  if (angleCount!=nAngles) {printf("Error: angleCount=%d!=%d=nAngles.\n",angleCount,nAngles); exit(0);}
}

/**********************************************************************************************
 * WriteLammpsForcefield -- writes "forcefield." file to be included into LAMMPS input script *
 **********************************************************************************************/
void WriteLammpsForcefield( FILE *fp ) 
{
  int i,j,n;
  real springRigidity_LAMMPS_units = springRigidity * J_IN_cal / Sqr(nm_IN_Angstrom);
  real angleRigidity_LAMMPS_units = angleRigidity * J_IN_cal;
  real radians_in_degrees = ( 180 / PI ); // conversion factor

  //  fprintf(fp, "pair_style hybrid lj/cut 12 lj/cut/coul/cut 12 dipole/cut 9\n\n");
  fprintf(fp, "pair_style dipole/cut 12\n");
  fprintf(fp, "bond_style harmonic\n");
  fprintf(fp, "angle_style cosine/squared\n");
  fprintf(fp, "\n");

  for ( n = 0; n < nTypes; n++ ) { // write masses
    fprintf(fp, "mass%5d%9.3f\n", n+1, mass[n]);
  }
  fprintf(fp, "\n");
  
  // pair_coeff for 'dipole/cut' style
  // syntax: pair_coeff iType jType ijEpsilon ijSigma cutoffLJ_optional cutoffElectro_optional
  for ( i = 0; i < nTypes; i++ ) { 
    for ( j = 0; j < nTypes; j++ ) {
      if ( i <= j ) { // LAMMPS restriction to avoid redundancy
	//      fprintf(fp, "pair_coeff%5d%5d\tdipole/cut%7.3f%7.3f\n", i+1, j+1, eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom); // hybrid style
	fprintf(fp, "pair_coeff%5d%5d\t%7.3f%7.3f\n", i+1, j+1, eps[i][j]*J_IN_cal, sig[i][j]*nm_IN_Angstrom);
      }    
    }
  }
  fprintf(fp, "\n");

  if (nSites!=nWaters) { // to prevent segmentation fault when simulating pure water systems
    // bond_coeff for 'harmonic' style
    // syntax: bond_coeff bondType rigidityConstant referenceLength
    // NOTE: in LAMMPS, this potential is defined as K(r-r0)^2, so K_LAMMPS = k_brahms/2
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHOL_PHOS, springRigidity_LAMMPS_units/2, 0.9*sig[CHOLINE_TYPE][PHOSPHATE_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", PHOS_GLY, springRigidity_LAMMPS_units/2, 0.9*sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLY_EST, springRigidity_LAMMPS_units/2, 0.9*sig[GLYCEROL_TYPE][ESTER_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", EST_CHAIN, springRigidity_LAMMPS_units/2, 0.9*sig[ESTER_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", CHAIN_CHAIN, springRigidity_LAMMPS_units/2, 0.9*sig[TAIL_TYPE][TAIL_TYPE]*nm_IN_Angstrom);
    //  fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", AMINE_PHOS, springRigidity_LAMMPS_units/2, 0.9*sig[PHOSPHATE_TYPE][AMINE_TYPE]*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYPLUS_PHOS, springRigidity_LAMMPS_units/2, (sig[PHOSPHATE_TYPE][GLYCEROL_TYPE]-sigma[GLYCEROL_TYPE]/4)*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTPLUS_CHAIN, springRigidity_LAMMPS_units/2, (sig[ESTER_TYPE][TAIL_TYPE]-sigma[ESTER_TYPE]/4)*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", GLYCHARGE_COM, springRigidity_LAMMPS_units/2, sigma[GLYCEROL_TYPE]/4*nm_IN_Angstrom);
    fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", ESTCHARGE_COM, springRigidity_LAMMPS_units/2, sigma[ESTER_TYPE]/4*nm_IN_Angstrom);
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 10, springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 11, springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 12, springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 13, springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
    /*   fprintf(fp, "bond_coeff%5d%7.3f%7.3f\n", 14, springRigidity_LAMMPS_units/2, 0.9*sig[][]*nm_IN_Angstrom); */
    fprintf(fp, "\n");

    // angle_coeff for 'cosine/squared' style
    // syntax: angle_coeff angleType rigidityConstant referenceAngle
    // NOTE: in LAMMPS, this potential is defined as K[cosTheta-cosTheta0]^2, so K_LAMMPS = k_brahms/2
    fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _phos_, angleRigidity_LAMMPS_units/2, refAngle_CholPhosGly*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _gly_, angleRigidity_LAMMPS_units/2, refAngle_PhosGlyEst*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piSat_, angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _cisUnsat_, angleRigidity_LAMMPS_units/2, refAngle_cisUnsat*radians_in_degrees);
    fprintf(fp, "angle_coeff%5d%8.4f%8.2f\n", _piDipCharges_, angleRigidity_LAMMPS_units/2, PI*radians_in_degrees);
    fprintf(fp, "\n");
  }
}
